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Polynomial approximation problems represent a class of specially structured problems which are frequently 
encountered in empirical curve-fitting. Two generators for creating such problems have been developed, 
implemented and used in the testing of discrete L\ approximation codes. Both generators permit automatic 
generation of problems with specified characteristics and (for one generator) having known, unique and 
controllable solutions. 

Key words: Algorithm testing; approximation; computational experiment; least absolute deviation; polynomial 
approximation; test problems 

1. Introduction 

Recent years have seen increased interest in least absolute deviation (L r ) data fitting, either as an 
alternative to or in conjunction with the usual least squares approach [11, L9]. 1 Estimation in theLj norm has 
certain desirable statistical properties (such as "robustness" [14, 16] when the underlying error distribution is 
long-tailed) and it can be carried out using reasonably efficient mathematical programming procedures. 

While a number of algorithms for discrete Lj approximation have recently been advanced [1, 3, 5, 6, 21, 
22], computational comparisons of specific implementations ("codes") for these approaches have not only 
been limited [3, 6, 22] but have often produced conflicting evidence. Resolution of these conflicts awaits the 
establishment of reliable and comprehensive methodology for testing such mathematical software. 

As a first step in developing a sound evaluation methodology for comparing mathematical programming 
codes, four L\ codes (representing a range of solution techniques and implementation strategies) were 
evaluated on particular classes of test problems [13]. The first group of test problems were "hand-picked" 
problems [20] — those created with a specified structure in mind or arising from actual applications. The 
second group consisted of pseudo-randomly generated problems, produced by means of a test problem 
generator [17] and representative of a range of general L x data fitting problems; such problems could be 
constructed with a number of controllable characteristics such as degeneracy, rank loss and optimal solution. 

The two generators P0LY1 and POLY2 described in this paper have subsequently been utilized to produce 
other types of problem classes for which L x approximation is appropriate. In particular, these are problems in 
which it is desired to obtain the best L x polynomial approximation to a specific function (POLY1) or to a set of 
discrete observations (POLY2). Polynomial approximation problems are frequently encountered in practice 
(e.g., among our hand-picked problems), and also enjoy an extensive theoretical basis [7, 18, 23]. Moreover, 
such problems are known to admit a range of "ill-conditioning" and numerical difficulties [10], notably 
evidenced in our previous testing efforts using 27 hand-picked problems [13] and in recent results obtained 
using the first of these two generators [8]. 

The use of generated test problems in code evaluation offers several advantages. Generators provide a 
virtually inexhaustible source of test problems, and a source that can be made relatively transportable from 
one computer to another. Since generators are able to provide test problems with controllable characteristics 
(such as size and structure), fairly well-defined classes of problems can be efficiently created. Moreover, 
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through the invocation of pseudo-random number generators, test problems can be "randomly" selected from 
such problem classes and valid statistical conclusions can be drawn about code performance with respect to 
these classes. Another desirable property of test problem generators, and one that is useful in judging the 
accuracy or correctness of a code, involves the ability to specify in advance the solution to a generated 
problem. 

Two generators are presented here for obtaining polynomial approximation test problems for L x curve- 
fitting. Both generators permit automatic generation of test problems with stipulated characteristics, and one 
of them (POLY2) allows the specification of known and unique solutions to the generated problems. The 
theoretical design of both generators is discussed in section 2, and the following two sections describe 
computer implementations of the theoretical design. The Appendices contain FORTRAN listings for both 
generators, written to be machine-independent. 2 



2. Design of the Generators 

The discrete linear L t approximation problem can be formulated as follows. Given a set of n observations on 
a single dependent variable y and each of m + 1 independent variables z , . . . , z m , find parameters /3 , 
. . . , j8 m that minimize 

n m 

i=\ 3=0 

The L i polynomial approximation problem , with which we will be concerned, is a restricted case of the above 
formulation: namely, minimize by choice of /3 = (/3o, . . . , fi m ) 

n m 

In the above, y { refers to the i-th observation on variable y and x { refers to the i-th observation on (the single 
independent) variable x. Equivalently, it is required to find a polynomial of degree m which best fits the given 
data, in the sense of minimizing the sum of absolute values of the residuals 

m 

?i = Ji ~2 PjXi- 
3=0 

A vector (3* which yields the minimum value </>* of (/>(/3) is termed a solution vector, with optimum objective 
function value (/)*. Unlike the case of L 1 polynomial approximation over a continuous interval [23, p. 38], the 
optimum solution vector /3* to a discrete Li polynomial approximation problem need not be unique. It is 
known [2, 4, 12] that such a solution vector can always be found which interpolates at least m + 1 of the data 
points; that is, at least m + 1 residuals e t are identically zero at the solution. A given problem exhibits 
degeneracy if more than ra + 1 residuals equal zero at the optimum. 

A number of algorithms for solving discrete L\ approximation problems make essential use of the above 
"interpolation" or "extreme point" result. Namely, these algorithms examine only basic solutions, where 
precisely ra + 1 residuals are zero, and move in a systematic way from one basic solution to another. Any 
such basic solution is defined by a set of ra + 1 indices (or rows) i e {l, 2, . . . , n} where e t = 0. For 
convenience, these indices are said to correspond to active constraints, while the remaining indices correspond 
to inactive constraints. Thus, in the absence of degeneracy the set N — {l, 2, . . . , n} can for any basic 
solution be partitioned as /V = /V° U N + U /V~, where /V° indicates the active constraints, /V + the inactive 
constraints with positive residuals, and /V - the inactive constraints with negative residuals. 



2 All FORTRAN programs have been checked for portability using the PFORT verifier (B. G. Ryder, Software-Practice and Experience, Vol. 4, 1974, 359- 
377) and conform to a subset of ANSI FORTRAN (X3.9-1966, Amer. Nat. Stand. Inst., New York, 1966). 
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The two generators for polynomial approximation problems discussed here will, for specified n and ra, 
produce the following test problem data: 

A = (*/), i = l,...,n; ; = m? (1) 

and 

y = (y f ), i=] ,». (2) 

The basic difference between the generators resides in the type of problem structure they simulate. 

The first generator, POLY1, is intended to model the situation where one wishes to approximate, over a 
discrete set, a continuous function f(x) by a polynomial of specified degree. Accordingly, the values y, in (2) 
are given by y, = /(#,-). Controllable features of this generator include: the function f{x), the real interval / 
over which f(x) is to be approximated, the number and distribution of observation points Xj e /, and the 
degree m of the approximating polynomial. A more detailed description of this generator is provided in 
section 3. It should be noted that POLY1 does not produce a known solution to the test problem created; 
therefore, other means are necessary in order to verify whether or not a given Li code has actually "solved" a 
generated problem. For this reason, we have also devised a computer routine to check optimality of code- 
produced solutions by examination of the Kuhn-Tucker conditions for an associated linear program [4, 17]. 

The second generator, POLY2, models the situation where an L x fit is required to a given discrete set of 
data (ji, *j), i = I, . . . , n. Unlike the case for P0LY1, the y, are not assumed to be generated via some 
known functional relation. Controllable features of this second generator include: the real interval of 
approximation /, the number and distribution of observation points x t e /, the degree of the approximating 
polynomial, and the statistical distribution of the residuals e r Moreover, one is also able to specify in advance 
a solution vector /3* which will be guaranteed to be the unique Lj solution to the generated problem. 

The approach used for generating a data set (A, y) with known optima] solution /3* will now be outlined. 
The matrix X in (1) is partitioned as 



X = 



A' " 

A + 

x- 



corresponding to indices i e /V° , /V + , and N~ respectively. Note that A° is a square (m + 1) X (m + 1) 
nonsingular matrix, with the proviso that the^/s fori e/V° are distinct. By a result proved in [17], a sufficient 
condition for the optimality and uniqueness of /3* in (A, y) is that 

2j x{ = 2*i x{, j = 0, . . . , m, (3) 

ieN + ieN- 

and that 

y = A/3* + e, (4a) 

where 

e t = if UN\ 

ei >0 if ieN + , (4b) 

d < if ieN~. 
It will now be shown how values x t can be found that satisfy (3). 



' For the first generator POLYl, one may also specify whether the polynomial has a constant term (j— m) or not (j= 1, . . . , m). 
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First, we note that fory = 0, equation (3) requires |/V + | = |W~|. For notational convenience, let r = |/V + | 
= |/V~| and redefine the x-s (i eN + ) as v u . . . , v r ; similarly, let the x-s (i eN~) be denoted as^!, . . . , 
w ,.. Also, write 

[v u . . . , v r \ = [iv u . . . , w r ] g (5) 



Zs Vp = 2j w p-> f° r h= 0,1, . . . , 6/. 

Then, following the development in [15], it is straightforward to verify that for any d 

[vi, . . . , v r ] Q = \wi, . . . , w r ] g => [vi, . . . , v r , w\ + d, . . . . w r + G^]^! 

(6) 

— [Wi, . . . , W r , V\ + d, . . . , V r + G?] g +i. 

Relation (6) can be used recursively to produce values x ( that satisfy (3). For example, since 

Mo = [l]o, 
equation (6) gives 

[0. 3]! = [1, 2]!, using d = 2, 
[0, 3, 5, 6] 2 = [1, 2, 4, 7] 2 , using d = 4, 

and so forth. Accordingly, the consecutive integers 0, 1, 2, . . . , 2 y+1 — 1 can be divided into two disjoint 
subsets {i>i, . . . , v r } and {10 15 . . . , w r }, with r = 2 Q , such that \p u . . . , v r ] a = [t^i, . . . , w r ] Q . The first 
subset consists of all integers I, < I <2 Q+1 — \ with an even number of ones in their binary expansion, and 
the second subset consists of all such integers with an odd number of ones. Dividing each v p and w p by 2 Q+l 
— 1 produces a set of numbers equally spaced on [0, 1], which in turn can be scaled and translated to yield 
a set of numbers {v p , w p : p — 1, . . . , r} equally spaced on an arbitrary interval / = [a, b]. Moreover, 
relation (5) still holds for the scaled and translated v p , w p . Finally, these values can be used in (3) as the 
required x-s for/ eN + and i eN~, respectively. 

More generally the recursive process based on (6) can be initiated, using any positive integer d , by means 
of 

[0]o = Mo 

and propagated using any d u d 2 , . . . , d Q to produce 2 Q+1 numbers satisfying (5), and therefore (3). We have 
chosen (see sect. 4) to generate the d-s from a uniform probability distribution. The resulting (scaled and 
translated) values for*, are no longer equally-spaced on /, but tend to be approximately normally distributed 
within the (finite) interval/. With probability one, all such Xj values generated are distinct. 

The above procedures form the theoretical basis for implementation of the second generator. Either 
equidistant x { or nonequidistant x f (corresponding to inactive constraints) can thus be generated. To ensure 
that (3) holds fory = 0, . . . , m, the number of observations n is required to be of the form 

n = m + 1 4- 2 m+k+ \ (7) 

where k is an integer, k > 0. Further details of P0LY2 are described in section 4. 
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3. Computer Implementation of POLY1 

This section describes a computerized version, in FORTRAN, for the first generator (POLY1) discussed in 
section 2. We will indicate those problem design characteristics available as user-specified input options to 
P0LY1, and then discuss how these inputs are utilized in creating test problem data. 

A variety of characteristics are available to the user as controllable options and are specified through input 
parameters to P0LY1. These options include: 

1. The number of observations n. 

2. The degree m of the approximating polynomial. 

3. The specific function /(jc) to be approximated. Ten such functions, 4 representing a variety of shapes 
likely to be encountered in practice, have been incorporated within the program (see table 1 and fig. 

i). 

4. The interval of approximation / = [a, 6]. 

5. A switch ICOL to indicate whether the approximating polynomial includes a constant term or not. 
(The latter case is useful in forcing the approximating polynomial to pass through the origin.) 

6. The location of the n observation points x x within /: either drawn from a uniform probability 
distribution or equally-spaced over/. 

7. A perturbation factor EPS, possibly zero, used in perturbing equally-spaced observation points*,-. 

8. Two initial starting seeds, utilized internally by a pseudo-random number generator. 

TABLE 1. Functions and Their Intervals of Approximation . 
fix) I = [a, b] 



1. e- x 8in* [0, 4] 

2. e*sin* [0, 4] 

3. e** n * [0, 7] 

4. e 2 */2x [.05, 1| 

5. 75*/[l + (7.5*) 2 ] [0, 2] 

6. \0xe~ ■** [0, 4] 

7. 1/[1 + (* - 2.5) 4 ] [0, 5] 

8. * lj3 [-1, 1] 

9. * lj2 [0, lj 

10. 1/[1 +x 4 ] [-2.5, 2.5J 

Given these input parameters, generation of the data sets can begin. The first step involves selecting the n 
observation points x t from the desired interval, according to the distributional form specified. The interval / 
= [a, b] is determined by the user (if CSWTCH =£ 0) through input parameters ENDL = a and ENDR = b, 
or (if CSWTCH = 0) by means of the default interval settings shown in table 1. The observation points are 
chosen according to a uniform probability distribution over / (if MSWTCH 4 1 0) or equally-spaced over / (if 
MSWTCH = 0). In the latter case, the user has the ability of perturbing the equally-spaced observation 
points #,•: namely, if/ = [a, b] then 

3d -> Xi + EPS-A-a,, i = 2, . . . , n - 1, 

where A = is the common subinterval length for equal spacing, the a ( are independent uniform 

n — 1 

random variates over [— 1, 1], and EPS is a user-supplied input representing the (maximum) fraction of the 
subinterval length A used for perturbation. If EPS = 0.0, no perturbation is performed; otherwise, a 
reasonable range for EPS is < EPS < V2. 

Next, the values y x = f(xi) are calculated, as well as the successive powers x{. If the user has specified a 
constant term in the approximating polynomial (ICOL = 1), then the above indexy has rangey = 0, 1, . . . , 
m. Otherwise (ICOL £ 1), a constant term is not included and j = 1, . . . , m. This stage of calculating 
successive powers xj is most sensitive to finite-precision arithmetic and to accumulated round-off error. 
Accordingly, such quantities are calculated in double-precision, using the recursion xj +l = Xi(x{). 

4 It will be noted that function 10 is simply a translation of function 7. It has been included as a separate entity to illustrate the effect of "scaling": namely, 
the successive powers .v,' produce more ill-conditioning when | x,| is larger (function 7) than smaller (function 10). 
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Finally, the problem data are returned in a single matrix PRBMAT. The first column of PRBMAT contains 
the vector y, and the remaining columns contain the design matrixX. In addition, the number of rows (NOBS) 
and the number of columns (NPAR) in matrix PRBMAT are returned to the calling routine, as well as a 
single-precision vector containing the observation points x { . The final seeds available from the random 
number generator are also returned in ISEED and JSEED, thus allowing the creation of different test problems 
(having the same input specifications) in successive calls of P0LY1 from a main routine. 

All error messages are written out using unit IOUT = 6. Changing this specific value of IOUT once at the 
beginning of the program permits other print unit numbers to be accommodated. Appendix A gives a complete 
listing of P0LY1, together with the subroutine UNIRAN [9], a well-tested uniform random number generator 
which it calls. 
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FIGURE 1. Graphs of functions. 
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Figure 2 illustrates the type of test problem that can be created by POLY1. The // = 50 data points (y h x t ) 
are derived from function 7 of table 1 over/ = [0, 5] and are shown here by the symbol "x"; in this case, the 
x-s were chosen to be equally-spaced (without perturbation) over /. By way of reference, the solid curve 
indicates the best L x fitting polynomial of degree m — 5 to these generated data. 

4. Computer Implementation of POLY2 

This section describes a computer implementation of the second generator (POLY2) for polynomial 
approximation problems. As indicated in section 2, an important feature of this generator is its ability to 
create data sets whose optimal Li solution vector can be specified in advance. Consequently, the correctness 
and accuracy of solutions produced by various Lj codes can be readily assessed using the known solutions to 
such problems. The fact that these solutions are guaranteed to be unique is also advantageous, since in the 
case of nonunique solutions it becomes difficult to compare the efficiency and accuracy of codes that reach 
correct, but different, solutions. 

A number of problem characteristics can be controlled by the user through appropriate input specifications. 
Options available in POLY2 include: 

1. The degree m of the approximating polynomial. 

2. The value k in (7), which together with in determines the number of observations n . 

3. The (unique) solution vector BETA to the generated problem. 

4. The interval of approximation / = [a, 6]. 

5. The location of the observation points x f corresponding to active constraints: either drawn from a 
uniform probability distribution or equally-spaced over/. 

6. The location of the observation points x f corresponding to inactive constraints: either equally-spaced 
or unequally-spaced over/. 



Y = l . .Q/( l + ( x - 2. 5,r ) 
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Y = 1.0/(1 + x H ) 
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Figure 1. (cont.) 
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Figure 2. Sample test problem created by POLY I for function 7, m — 5, and n = 50. 
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7. A perturbation factor EPS, possibly zero, used in perturbing equally-spaced x t associated with active 
constraints. 

8. The distribution of the residuals e,- associated with inactive constraints. 

9. Two initial starting seeds, utilized internally by pseudo-random number generators. 

The first step in developing a test problem is to select the solution vector BETA, if the user has not chosen 
to input a vector to POLY2. The solution vector is randomly selected, in this case, from a uniform probability 
distribution over an interval specified by the user (BSW = 1) or defined within the program (BSW J= 1). 

Next, the observation points x ( are generated from the interval / according to the desired distribution. 
Interval / is either specified on input (CSW = 1) or assigned a default setting within the program (CSW =/= 1). 
The x-s corresponding to active constraints are derived from a uniform probability distribution over /(MSW 
= 1) or are equally-spaced over /(MSW f 1). In this latter case, the x { can be perturbed using the factor 
EPS, exactly as described for POLY1. The observation points x t corresponding to inactive constraints are 
generated so that equation (3) is satisfied. Such x f are either equally-spaced (ISW = 1) or unequally-spaced 
(ISW £l) over the interval /; see section 2 for details of how these two methods of generation are achieved. 

Also, the entries xj of the design matrix are successively calculated fory = 0, . . . , m using appropriate 
care to preserve numerical accuracy. At the same time, the residuals e h i e /V + or i e/V~, are generated from 
either a uniform (RSW = 1) or a normal distribution (RSW ^l), are given the appropriate sign, and are then 
used to calculate y ( from (4a) — (4b). The optimum objective function value SUMRES is calculated in double- 
precision, based on these residuals. The input parameter RESID indirectly controls the magnitude of this 
objective function value. Namely, RESID acts as a scale factor with respect to the size of the residuals, which 
are either selected from a uniform distribution on [ — RESID, RESID], whose standard deviation is thus 
RESID/^3, or from a normal distribution having mean and standard deviation RESID. 

Finally, the problem data are returned as before in a single matrix PRBMAT. The first column of PRBMAT 
contains the y vector, and the remaining columns contain the design matrix X. The number of rows (KKM1) 
and the number of columns (M2) are also returned, as well as BETA, SUMRES and LOC (a vector containing 
the m + 1 row locations of the active constraints). The final seeds available from the last call to the random 
number generator are returned in ISEED and JSEED. 
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All error messages are written out using unit IOUT, which can be easily modified to accommodate differing 
print unit assignments. Appendix B gives a complete listing of P0LY2, together with the subroutines it calls: 
SORTM, SET1, SET2, EQUAL, NORRAN, SORTP and UNIRAN. The call structure of these subroutines is 
shown in figure 3; the latter three subroutines are fairly well-tested and are described further in [9]. 

Figures 4 and 5 illustrate the types of data sets that can be produced using POLY2. Both show as dots 134 
data points ( y i . v,) together with the best L x fitting polynomial of degree 5 (solid curve). In figure 4, the % t are 
generated via the equally-spaced option (ISW = 1), while in figure 5 they are generated via the unequally- 
spaced option (ISW T 1). 




Figure 3. Subroutine call structure for P0LY2. 




X 
FIGURE 4. Sample test problem created by POLY2: n = 134, m = 5, equally-spaced %\. 
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X 

FIGURE 5. Sample test problem created by POLY2: n = 134, m = 5, unequally-spaced x { . 
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6. Appendix A. 

Listing of FORTRAN Subroutines for First Generator: 

POLY1 
UNIRAN 



SUBROUTINE PC LY 1 ( NFUNC * M .N . PkBM AT • CSWTCH. MSWTCH ,ENDL • 
+ ENDR t I COL 9 NOB S . NP A~ . I SEED • J SEED ♦ X . N3: W . E PS ) 
C 

Q *********************************** ******** *************************** 

c 

c description: 

c 

C THIS SU3R0UTINE GENERATES DISCRETE POLYNOMIAL APP^CX IM AT I ON 

C PROBLEMS WHICH CAN BE USEO IN THE TESTING AND EVALUATION OF 

C ALGORITHMS FOR LI (LEAST ABSOLUTE DEVIATION) CUFVE FITTING* 

C THE USER CAN SPECIFY THE FOLLOWING CHARACTERISTICS OF THE 

C GENERATED DATA SETS: 

C 

C 

C * PROBLEM SIZE 

C * THE SPECIFIC FUNCTION TO BE APPROXIMATED 

C * THE REAL INTERVAL OVER WHICH THE FUNCTION IS APPROXIMATED 

C * THE DISTRIBUTION OF THE OBSERVATION POINTS IN THAT INTERVAL 

C 

C 

C THE POLYNOMIAL APPROXIMATION PROBLEM CONSIDERED HEF.E IS TO 

C APPROXIMATE (IN THE LI NORM) A GIVEN FUNCTION BY A POLYNOMIAL 

C OF SPECIFIED DEGREE. OVER A DISCRETE SET. F03 1=1. N X(IJ IS 

C A SELECTED POINT FROM SOME INTERVAL AND Y{I) IS THE CQf RES- 

C PONDING FUNCTIONAL VALUE AT X(I). IN THE PROBLEMS GENERATED 

C BY THIS SUBROUTINE. THE CONSTANT TERM OF THE POLYNOMIAL CAN BE 

C INCLUDED OR EXCLUDED (IN THE L.ATTER CASE FORCING THE APPRCX- 

C IMATING POLYNOMIAL TO PASS THROUGH THE ORIGIN). 

C 

C FOR CONVENIENCE THE Y VECTOR AND THE DESIGN MATRIX. WHICH 

C CONTAINS SUCCESSIVE POWERS OF THE OBSERVATION POINTS X ( I .) . 

C ARE BOTH RETURNED IN A SINGLE MATRIX •PRBMAT'. THE FIRST 

C COLUMN OF PRBMAT CONTAINS THE Y VALUES AND THE REMAINING 

C COLUMNS CONTAIN THE DESIGN MATRIX© 

C 

c note: 
c 

C * THE X(I) VALUES ARE AR&ANGED TO BE IN INCREASING ORDER. 

C FOR 1 = 1. N. IN PARTICULAR. X(l) IS ALWAYS THE LEFT-HAND 

C ENDPOINT OF THE INTERVAL. AND X(N) IS ALWAYS THE R IGHT- 

C HAND ENDPOINT. 

C 

c note: 
c 

C * SINCE AN OPTIMAL SET OF COEFFICIENTS FOR THE APPROXIMAT- 

C ING POLYNOMIAL IS NOT PRODUCED BY THIS GENERATOR. VERIFI- 

C CATION OF AN OPTIMAL SOLUTION SHOULD BE MADE BY SOME MEANS 

C (E.G.. EXAMINING THE KUHN-TUCKER CONDITIONS). 

C 

C*** *********************************************** ********** ********** 

c 

c input arguments: 
c 

C NFUNC - DETERMINES THE SPECIFIC FUNCTION BEING APPROXIMATED 

c M DEGREE OF THE APPROXIMATING POLYNOMIAL 

C N NUMBER OF OBSERVATION POINTS X(I) 
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C CShTCH- INPUT SWITCH USED F OR INTERVAL SELECTION, IF 

C = 1 USES ENDPOINTS SUPPLIED BY THE USER TO DEFINE 

C THE INTERVAL 

C =0 SELECTS THE INTERVAL USING DEFAULT SETTINGS 

C MSKTCH- INPUT SWITCH FOi- THE DISTRIBUTION GF THE OBSERVATION 

C POINTS. IF 

C = 1 PRODUCES JOINTS FFOM A UNIFORM ^F.GbABILITY DIS- 

c TfiiBUTiCN ove; THE INTEhVA^. 

C = C PRODUCES POINTS EQUALLY SPACED OVER THE INTERVAL 

C ( THOUGH POSSIBLY PE^TU r BED BY A RANDOM AMOUNT 

C If EPS > ) 

C ENOL — USER SPECIFIED LEFT ENDPOINT OF THE INTERVAL 

C ENQR USER SPECIFIED RIGHT ENDFCINT CF THE INTERVAL 

C I COL INPUT SWITCH. IF 

C = 1 P.'OuUCES AN Li PKCJLEM WITH A CONSTANT COLUMN 

C =0 PRGDUCES AN LI PROBLEM W/C A CONSTANT COLUMN 

C ISEED - Fl'ST OP T*C INITIAL SEEDS FOR THE aANDOM NUMBER 

C GENERATOR 

C JSEED - SECOND OF T wO INITIAL SEEDS FOE THE RANDOM NUMBER 

C GENERATOR 

C N&OW — MAXIMUM NUMBL.K OP CBSEfVATION POINTS 

C EPS FRACTION OF THE SUI3INTERVAL LENGTH USED IN PER- 

C TUBBING EQUALLY-SPACED OBSERVATION POINTS. IF 

C MSWTCH=0 

C 

c note: 
c 

C * WHEN EPS = 0. NO PtRTURBATICN IS PERFORMED. OTHERWISE. A 

C REASONABLE RANGE IS 0<E D S<0.5 • 

C 

C * THE DEFAULT VALUE FOP EACH OF THE SWITCHES I COL. CSWTCH. 

C MSWTCH IS ZERO. 

C 

c output arguments: 
c 

C PrdMAT- THE GENERATED PROBLEM MATRIX. WITH NOBS SO* AND NFAr 

C COLUMNS 

C NOBS — NUMBER OF ,\CWS IN PP.BMAT 

C NPAR — NUMBER OF COLUMNS IN PROMAT 

C ISEED - FISST OF TWO FINAL SEEDS AVAILABLE FROM THE RANDOM 

C NUMBER GENERATOR 

C JSEED - SECOND OF TWO FINAL SEEDS AVAILABLE FkCM THE RANDOM 

C NUMBER GENERATOR 

C X VECTOR Of OBSERVATION POINTS. OF DIMENSION NOBS 

C 

c 

c restrictions: 
c 

C N MUST BE GT 1 

C N MUST BE LE NlOW 

C M ^JST BE GT 

C N CANNOT BE LT M 

C NFUNC MUST BE BETWEEN 1 ANU 10, INCLUSIVE 

C 

(^*********** *************************¥********************************* 

c 

C A VAUETY OF FUNCTIONAL FC1MS AND SHAPES ARE CUi F ENTLY AVAILABLE: 

C 
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c 
c 

c 
c 


NFUNC FUNCTION DEFAULT INTERVAL 


1 EXP(-X)*SIN(X ) [J.0.4,0] 


c 


2 EXP{X)*S IN(X ) CO. 0,4.0] 


c 


3 EXP(X*SIN<X) ) CC. 0,7.0] 


c 


4 o5*£X?{2e*X)/X C. 05, 1.0] 


c 


5 75.*X/< 1 . + (7.5*X )**2) CO. 0,2.0] 


c 


6 10.*X*EXPC-.5*X) CO. 0,4. J] 


c 


7 l.0/(loJ+(X-2e5)**4) CO. 0,5.0] 


c 


8 CBRT (X) C -1 .3, 1.0 ] 


c 


9 SORTC X) C CO. 1.0 ] 


c 

c 
c 
c 
c 


10 U0/(U0HX**4>) C-2.5.2.5] 


SUBPROGRAMS USED! 


SIN FORTRAN LIBRARY FUNCTION 


c 


EXP FOf.T^AN LIBRARY FUNCTION 


c 


CBRT — FORTRAN LIBRARY OR USER- SUPPL IED FUNCTILN 


c 


v TO COMPUTE CUBE r COTS 


c 


SQRT FORTRAN LIBRARY FUNCTION 


c 


FLOAT - FOrJFAN LIB3A3Y FUNCTION 


c 


SORT — SUBROUTINE USED TO SORT THE ELEMENTS OF A VECTOR GF 


c 


DIMENSION N INTO ASCENDING ORDE^ 


c 


UNIRAN SUBROUTINE TO GENERATE PSEUDO-RANDOM NUMBERS FROM 


c 
c 
c 
c 
c 
c 
c 


A UNIFORM DISTRIBUTION Ov/EA [0,1 ] 


LANGUAGE: ANSI FC^T.'AN (1966) 


references: 


K.LoHOFFMAN & Do <o SH 1 E; r • ■ A TEST PROBLEM GENERATOR FO? 


c 


DISCRETE LINEAR LI APPROXIMATION PROBLEMS^, TO APPEAR IN 


c 


TRANS* ON MATH© SCFTfcA^E 


c 


PoDoDOMICH & Oo-oSHIEF , •GENEfATOrS FOR DISCRETE POLYNOMIAL 


c 
c 

c 
c 

c 


LI APPROXIMATION PROBLEMS* 


WRITTEN BY: 


PAUL D. DOMICH 


c 


CENTER FD~l APPLIED MATHEMATICS 


c 


NATIONAL BUREAU OF STANDARDS 


c 

c 

c 

c 

c** 

c 


WASHINGTON D0C0 2023^ 


JANUARY, 1979 


*********** ********************* *************************** ********** 


DOUBLE PRECISION UNI I NC , EPSA , DA , D^aOD 




REAL PRDMATCNROW, 1 ) .ENDSC 10,2) ,X(NRCW) 




INTEGER CSa/TCH 


f" 


LOGICAL CONST 


c 

C***** DEFAULT SETTINGS FOR THE INTERVAL SPECIFICATIONS 


c 


DATA ENDS( 1,1) • ENDSC 2, I ) ,ENDS( 3,1) , ENDSC 4, 1 ) , ENDSC 5 • 1 ) , 




+ ENDS (6 ,1 ) , ENDSC 7, 1 ) , ENDSC 3, 1 ), ENDSC 9, 1 ), ENDSC 1C ,1 ), 




+ ENDSC 1,2) ,END5C2.2) ,ENDS( 3* 2). ENDSC 4 .2) , ENDSC 5, 2) • 



468 



+ ENDS(£t2)tEN0S{7,2).ENDS(8,2) *ENDS(9.2 ) .ENDSC 10t2)/ 

L 0©0» OoO* Oou» o05. OoOt OoOt 0©0. -UO« OoO* -2©5t 

R 4.0, 4.0. 7.c. 1.0 . 2.0. 4.3. 5.0. l.Ci 1.0* 2.5/ 

C 1234567 8910= NFUNC 

C 

C***** CHECK INPUT DATA FOi\ CONSISTENCY 
C 

I 0UT = 6 

IF(N .LE. 1 .OR. N .LT. M) GO TC 2001 

IMN oEUe M ©AND© ICOL o E Q© i) GOTO 2001 

IF(M .LT. 1) GOTO 2C01 

ITiN ©GT© N'OW) GCTC 2002 

IF(NFUNC .LT. I .OK. NFUNC .GT. 10) GOTO 2003 

IF(NFUNC .EG. 4 .AND. ENOL .LE. 0.0 .AND. ENDR . GE . 0.0) GCTC 2004 

IF (NFUNC .EQo 9 • AND© { END7. oLTo OoO »OR« LNDL ©LT© 0.0) ) GCTO 200£ 
C 

C***** CALCULATION OF THE INTERVAL 3Y SPECIFICATION Ck DEFAULT 
C 

IF(CSWTCH oEQc 1) A=ENDL 

IF(CSWTCH .EQ. 1) E = ENDF 

IF(C5WTCH ©NE© 1) A=ENDS(NFUNC t 1 ) 

IF(CSWTCH .NE. i) B=ENDS ( NFUNC . 2 ) 

IF(A .LE. B) GOTO 4 

HGLD=A 

A = B 

B=HOLD 

4 XINT=3-A 
C 

C***** CALCULATE THE REQUIRED NUMBER OF COLUMNS IN THE DESIGN MATRIX 
C 

IF (ICOL .EQ. 1) NVAR=M+1 

IF (ICOL ©EQ. 1) CGNST = oTf.UEo 

IF (ICOL .NL. 1) NVAR=M 

IF (ICOL cNEc I) CuNST=oF ALSEo 
C 

C***** CALCULATE THE OBSERVATION POINTS WITHIN THE INTERVAL CA.B3 
C 

CALL UNI.IANCN© i »X. I SEED, J SEED) 

X( 1 )=A 

X(N)=B 

K = N- 1 

IF(K©EQ© 1 ) GCTO 8 

IF (MSWTCH .EQ. i) GOTO 6 
C 

C***** CALCULATE THE SUBINTERVAL LENGTH FOR EQUALLY-SPACED POINTS. 
C POSSIBLY MODIFIED BY A PERTURBATION FACTOR (EPSA). 

C 

UNI INC=XINT/FLQAT( N- 1 ) 

EPSA=EPS*UNI I NC 

DO 5 1=2. K 

X( I )= A+FLCAT ( 1-1 )*UNI INC + EPSA*(2©0*X( I )-l ©0) 

5 CONTINUE 
GOTO 5 

C 

C***** CALCULATE THE OBSERVATION POINTS USING A UNIFORM DISTRIBUTION 

C AND SG..T THE POINTS INTO ASCENDING OHDEP 

C 

6 DO 7 I=2.K 

X( I >=XINT*X( I ) +A 

7 CONTINUE 

CALL S0~T(X,N,X) 



469 



c***** 

c 



10 



20 
30 

40 

50 

60 

7C 

80 

90 

100 
1003 
C 

c 



CALCULATION OF THE Y VECTOR 

DO 1030 1 = 1. N 

A=X( I ) 

GOTO( 10.20.30.40,5 0.6 0.70.30.90. 100) # MFUNC 

PRBMATC I . 

GOTO 1003 

PRBMATC I . 



GOTO 1003 
PRBMATC I , 
GOTO 1003 
PRBMATC I . 
GOTO 1003 
PRBMATC I • 
GOTO 1003 
PR3MATCI . 
GOTO 1003 
PRBMATC I, 
GOTO 1003 
PRBMATC I • 
GOTO 1003 
PRBMATC I 9 
GOTO "1003 
PRBMATC I t 
CONTINUE 

CALCULAT 

DA=A 

IFC CONST 

IFC •NOT. 



)=EXPC~A)*SINC A) 

>=EXPCA**SINCA> 



)=LXPCA*SINC A) ) 



)=.5*EXP(2.0*A)/A 



)=75.0*A/C I .0+C 7.5* A )**2) 



>=10.0*A*EXPC-.5*A) 



)-l»O/Cl«0>( A-2.5 )**4) 



)=CBRTC A) 



i=SQRTC A) 



) = 1*0/C 1«0+C A**4) ) 



ON CF THE ENTRIES IN THE DESIGN MATRIX 



) OPRCO=U0 
CONST I CPROD=A 



1004 
1030 

C 

c 



PRBMATC I .2 )-OPROD 

IFCNVAR.EQ. 1 ) GOTO 1 330 

DO 1004 J=2.NVA^ 

DPRCD=DPRCD*DA 

FSBMATC I • J + l )=DPROD 

CONTINUE 

CONTINUE 

DEFINE ROW AND COLUMN DIMENSIONS OF THE PROBLEM MATRIX CPKBMAT) 



N.M 



OF N 
NSOW 



N0B5-N 
NPAR=NVAR+1 
GOTO 3000 
20 1 WRITEC I0UT.9C 1 ) 

901 FORMAT! IH0.21HTHE INPUT VALUES N 
+ 18H ARE NOT FEASIBLE) 

GOTO 3000 

2002 WRITEC IOUT. 902) N.NROW 

902 FORMATC 1H0. 17HTHE VALUE 
+ 21H DIMENSION OF 

GOTO 3000 

2003 WRITEC IOUT. 903) NFUNC 

903 FO&MATC 1H0.22HTHE VALUE 
GOTO 3000 

2004 WRITEC IOUT. 904)NFUNC 

904 FCFMATC 1H0.30H THE INTEi 
GOTO 300 

2005 WRITEC IOUT. 9C4INFUNC 
3000 CONTINUE 

RETURN 
END 



=.I5.10H 



AND M 



IS . 



=. I5.21H 
= . I 5 ) 



EXCEEDS THE MAXIMUM. 



OF NFUNC 



=.15,18H 



IS NOT IN rANGE) 



VAL SET FOR FUNCT I ON . I 3 . 1 1 H IS INVALID) 



470 



S UBROUT I NE UN I RAN ( N . I START .X.I SEED . JSEED ) 
C 

C PURPOSE THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF S I ZE N 

C FROM THE UNIFORM (RECTANGULAR) 

C DISTRIBUTION ON THE UNIT INTERVAL (0.1). 

C THIS DISTRIBUTION HAS MEAN = 0.5 

C AND STANDARD DEVIATION = SQRTd/12) = 0.28867513. 

C THIS DISTRIBUTION HAS THE PROBABILITY 

C DENSITY FUNCTION FIX) = 1. 

C INPUT ARGUMENTS N = THE DESIRED INTEGER NUMBER 

C OF RANDOM NUMBERS TO BE 

C GENERATED. 

C ISTART * AN INTEGER FLAG CODE WHICH 

C C IF SET TO 0) WILL START THE 

C GENERATOR OVER AND HENCE 

C PRODUCE THE SAME RANDOM SAMPLE 

C OVER AND OVER AGAIN 

C UPON SUCCESSIVE CALLS TO 

C THIS SUBROUTINE WITHIN A RUN; OR 

C (IF SET TO SOME INTEGER 

C VALUE NOT EQUAL TO • 

C LIKE. SAY. 1* WILL ALLOW 

C THE GENERATOR TO CONTINUE 

C FROM WHERE IT STOPPED 

C AND HENCE PRODUCE DIFFERENT 

C RANDOM SAMPLES UPON 

C SUCCESSIVE CALLS TO 

C THIS SUBROUTINE WITHIN A RUN. 

C OUTPUT ARGUMENTS — X = A SINGLE PRECISION VECTOR 

C (OF DIMENSION AT LEAST N) 

C INTO WHICH THE GENERATED 

C RANDOM SAMPLE WILL BE PLACED. 

C OUTPUT — A RANDCM SAMPLE OF SIZE N 

C FROM THE RECTANGULAR DISTRIBUTION ON (0.1). 

C PRINTING NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 

C RESTRICTIONS THERE IS NO RESTRICTION ON THE MAXIMUM VALUE 

C OF N FOR THIS SUBROUTINE. 

C OTHER DATAPAC SUBROUTINES NEEDED NONE. 

C FORTRAN LIBRARY SUBROUTINES NEEDED NONE. 

C MODE OF INTERNAL OPERATIONS SINGLE PRECISICN. 

C LANGUAGE ANSI FORTRAN. 

C COMMENT ******************** 

C THE PARAMETERS OF THE CONGRUENTI AL GENERATOR 

C USED HEREIN ARE OF NECESSITY STRONGLY MACHINE- 

C DEPENDENT IN TERMS OF THE QUALITY OF THE 

C RANDOM NUMBERS PRODUCED. KNOWN GOOD RESULTS 

C HAVE BEEN OBTAINED ON THE UNIVAC 1108 (36 BIT) 

C COMPUTER. KNOWN POOR RESULTS HAVE BEEN OBTAINED 

C ON THE CDC 3300 (48 BIT) COMPUTER. 

C IT IS STRONGLY RECOMMENDED THAT UPON 

C IMPLEMENTATION OF UNIRANt THE OUTPUT 

C FROM THIS SUBROUTINE SHOULD BE CHECKED FOR INDEPENDENCE 

C AND UNIFORMITY. DATAPAC SUBROUTINES PLCTX. 

C PLOTXX. PLOTU. RUNS. TIMESE. HIST. TAIL. LOC. AND 

C SCALE MAY BE USEFULLY EMPLOYED TO DO SUCH 

C CHECKING. SUCH CHECKING IS ESPECIALLY IMPORTANT 

C IN LIGHT OF THE FACT THAT OTHER DATAPAC RANDOM 

C NUMBER GENERATOR SUBROUTINES (NORRAN — NORMAL. 

C LOGRAN — LOGISTIC. ETC.) ALL MAKE USE OF INTERMEDIATE 

C OUTPUT FROM UNIRAN. 
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c ******************** 

C REFERENCES KRONMAL. 'EVALUATION OF THE PSEUDO-RANDOM 

C NORMAL NUMBER GENERATOR 1 • JOURNAL OF THC 

C ASSOCIATION FOR COMPUTING MACHINERY * 1964, 

C PAGES 357-363. 

C HAMMERSLEY AND HANDSC0M8. MONTE CARLO METHODS. 

C 1964, PAGE 36 • 

C — JOHNSON AND KOTZ. CONTINUOUS UNIVARIATE 

C DISTRIBUTIONS 2. 1970, PAGES 57-74* 

C WRITTEN BY — JAMES J. FILLIBEN 

C STATISTICAL ENGINEERING LABORATORY <205.03) 

C NATIONAL BUREAU OF STANDARDS 

C WASHINGTON, D. C. 20234 

C PHONE: 301-921-2315 

C ORIGINAL VERSION — JUNE 1972. 

C UPDATED AUGUST 1974. 

C UPDATED SEPTEMBER 1975. 

C UPDATED --NOVEMBER 1975. 
C 

C 

DIMENSION X(l ) 
C 

DATA DIV.K1 .K2.K3. I 1,1 2/34359 738368. 0,4097 .129. 221 82922 491. 208673 5 
10019.18575103187/ 
C 

IPR=6 
C 

C CHECK THE INPUT ARGUMENTS FOR ERRORS 
C 

IF<N.LT.l IGOT050 

G0T090 
50 WRITE<IPR. 5) 

WRITE< IPR.47IN 

RETURN 
9 CONTINUE 
5 FGRMATUH . 91H***** FATAL ERROR — THE FIRST INPUT ARGUMENT TO THE 
1 UNIRAN SUBROUTINE IS NON-POSITIVE *****) 
47 FORMATIIH . 35H***** THE VALUE OF THE ARGUMENT IS .18 .6H *****) 

C 

C- START POINT — 

C 

I1=IS£ED 

I 2=JSEED 

IF( ISTART.NE.0)GOTO150 

11=20867350019 

12=18575103187 
C 

150 D0100I=1.N.2 

I 1=IABS< I1*KH-1> 

I 2=IABSC <I2*K2+K3J*K2+K3) 

Ul=FLOAT<Il)/DIV 

U2=FL0AT( I2)/DIV 

XII *=U1 

IFC I.EQ.NJGOTOIOO 

XCI+1 |=U2 
100 CONTINUE 
C 

ISESD=U 

JSEED=I2 

RETURN 

END 
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7. Appendix B. 
Listing of FORTRAN Subroutines for Second Generator: 

POLY2 

SET1 

SET2 

UNIRAN* 

NORRAN 

EQUAL 

SORTM 

SORTP 

* For a listing of UNIRAN, refer to Appendix A. 
SUBROUTINE PCLY2(W,K,NRCW,JSW, BETA , BSW » BLEFT t BR [ GHT , CS W . CLEFT , 

+ CRIGHTt I5W,MSW,EPS.kSW,RESID, I SEED, J SEED , PkBMAT, SUMREStKKMl , M2, 
+ LOC ) 
C 

C ****** ************************ **************** ************************* 

c 

c description: 

c 

C THIS SUBROUTINE GENERATES DISCRETE PCLYNOMIAL APPROXIMATION 

C PROBLEMS WHICH CAN BE USED IN THE TESTING AND EVALUATION OF 

C ALGORITHMS FOR LI (LEAST ABSOLUTE DEVIATION) CUf.VE FITTING. 

C THE USER CAN SPECIFY THE FOLLOWING CHARACTERISTICS OF THE 

C GENERATED DATA SETS: 

C 

C * PROBLEM SIZE 

C * THE UNIQUE SOLUTION VECTOR TO THE GENERATED LI PROBLEM 

C * THE REAL INTERVAL CVLF WHICH THE APPROXIMATION IS DEFINED 

C * THE DISTRIBUTION OF THE OBSEKVATICN POINTS WITHIN THE 

C INTERVAL 

C * THE DISTRIBUTION Ub THE RESIDUALS 

C 

C THE POLYNOMIAL APPROXIMATION P.300_EM CONSIDERED HE ?E IS TO 

C APPROXIMATE (IN THE LI NORM) A DISCRETE SET OF DATA VALUES Y 

C BY A PLLYNCMIAL OF SPECIFIED DEGREE OVE.n SOME INTERVAL. FOR 

C I=1,KKM1 X(I) IS AN OBSERVATION POINT IN THE INTERVAL ANO 

C Y(I) IS THE DATA VALUE CORRESPONDING TO X(I). THE VALUES FOR 

C X(I) AND Yd) A?,E SELECTED IN SUCH A MANNED TO GUARANTEE CPT- 

C IMALITY OF THE SOLUTION VECTOR BETA IN THE GENERATED PROBLEM, 

C 

C FOR CONVENIENCE, THE Y-VALUES AND THE DESIGN MATRIX CONTAINING 

C THE APPROPRIATE POWERS OF THE OBSERVATION POINTS X(I) ARE BOTH 

C RETURNED IN A SINGLE MATr.IX •P^BMAT'© THE FI^,ST COLUMN OF P^BMAT 

C CONTAINS THE Y-VALUES AND THE REMAINING COLUMNS CONTAIN THE DESIGN 

C MATRIX. 

C 

c note: 
c 

C * A CONSTANT TERM IS INCLUDED IN THE APPROXIMATING POLYNOMIAL 

C 

C * THE OBSERVATION POINTS ARE ARRANGED IN INCREASING ORDER, BY 

C ROW, WITHIN THE GENERATED PROBLEM MATRIX 

C 

C * THE PROCEDURE GIVEN HERE GENERATES AN LI APPROXIMATION 

C PROBLEM HAVING A KNOWN AND UNIQUE SOLUTION VECTOR BETA. 

C TWO METHODS OF GUARANTEEING THESE PFCPExTIES ARE PROVIDED 

C IN SUBROUTINES SETl AND SET2 WHICH GENERATE THE OBSERVA- 

C TION POINTS FCS THE INACTIVE CONSTRAINTS. THE METHOD USED 

C IN SUBROUTINE SET2 UTILIZES A RANDOM NUMBER GENERATOR AND 

C THU^» CAN CREATE DIFFERENT LI PROBLEMS IN SUCCESSIVE RUNS 

C HAVING THE SAME INPUT SETTINGS, WHEFEAS THE SAME PROBLEM 

C WILL BE REPRODUCED WHEN USING SUBROUTINE SETl. 

473 



c 

C******************************************* **** ***************** ****** 

c 

c input arguments: 
c 

C m DEGREE OF THE APPROXIMATING POLYNOMIAL 

C K INTEGER USED IN DETERMINING THE NUMBER OF INACTIVE 

C CONSTRAINTS TO THE PPGBLEM 

C NROW — MAXIMUM NUMBER OF OBSERVATION POINTS 

C JSW IN=>UT SWITCH USED IN THE SELECTION OF THE SOLUTION 

C VECTOR BETA, IF 

C = I PRODUCES A SOLUTION VECTOR WITH ELEMENTS DRAWN 

C FROM A UNIFORM DISTRIBUTION 

C - USES THE SOLUTION VECTOf SPECIFIED EXPLICITLY BY 

C THE USER 

C BETA USER SPECIFIED SOLUTION VECTOR TO THE GENERATED PROB- 

C LEM (IF JSw=0> OF DIMENSION M+l 

C BSW INPUT SWITCH USED TO DEFINE THE INTERVAL FROM WHICH 

C ELEMENTS OF THE BETA VECTOR ARE SELECTED (WHEN JSW=1) 

C IF 

C = I USES ENDPOINTS SUPPLIED BY USER TO DEFINE THE 

C INTERVAL 

C =0 SELECTS ENDPOINTS BY DEFAULT SETTINGS 

C BLEFT - LEFT ENDPCINT OF THE INTERVAL USED IN SELECTING 

C ELEMENTS OF THE SOLUTION VECTOR BETA (IF JSW=3SW=1) 

C BRIGHT- FIGHT ENDPOINT OF THE INTERVAL USED IN SELECTING 

C ELEMENTS OF THE SOLUTION VECTOF: BETA (IF JSW~BSW=1) 

C CSW INPUT SWITCH USED TO DEFINE THE INTERVAL FROM WHICH 

C OBSERVATION POINTS ARE SELECTED, IF 

C = 1 USES ENDPOINTS SUPPLIED BY USE? TO DEFINE THE 

C INTERVAL 

C - SELECTS ENDPOINTS BY DEFAULT SETTINGS 

C CLEFT - LEFT ENDPOINT OF THE INTERVAL USED IN SELECTING US- 

C SE3VATI0N POINTS (IF CSW=1) 

C CRIGHT- RIGHT ENDPCINT CF THE INTERVAL USED IN SELECTING OB- 

C SERVATION POINTS (IF CSW=1) 

C ISW INPUT SWITCH WHICH DEFINES THE METHOD OF GENERATING 

C OBSERVATION POINTS FOA THE INACTIVE CONSTRAINTS, IF 

C =1 PRODUCES EQUIDISTANT OBSERVATION POINTS BY 

C SET1 

C =0 PRODUCES NON-EQUIDISTANT OBSERVATION POINTS BY 

C SET2 

C MSW INPUT SWITCH WHICH SPECTFIFS DISTRIBUTION OF 

C OBSERVATION POINTS TO,; THL ACTIVE CONSTRAINTS* IF 

C = 1 PRODUCES OBSERVATION POINTS DRAWN FROM A 

C UNIFORM DISTRIBUTION 

C =0 PRODUCES OBSERVATION POINTS EQUALLY SPACED 

C CVE1 THE INTERVAL (THOUGH POSSIBLY PEFTUftBED 

C BY A RANDOM AMOUNT IF EPS > 0) 

C EPS FRACTION CF THE SUBINTERVAL LENGTH USED IN PEF,- 

C TURBING EQUALLY SPACED OBSERVATION POINTS FOR 

C THE ACTIVE CONSTRAINTS (IF MSW=0) 

C RSW INPUT SWITCH WHICH SPECIFIES THE DISTRIBUTION 

C OF THE (NONZERO) RESIDUALS ASSOCIATED WITH THE IN- 

C ACTIVE CONSTRAINTS, IF 

C =1 PRODUCES RESIDUALS DRAWN FrvCM A UNIFORM DIS- 

C TP I OUT I ON 

C =0 PRODUCES RESIDUALS DRAWN FROM A NORMAL DIS- 
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C TRIBUTIGN 

C ?*ESID - STANDARD DEVIATION OF THE RESIDUAL DISTRIBUTION 

C (WHEN RSW=0>. CR ABSOLUTE VALUE OF THE UPPER 

C AND LOWER LIMITS OF THE UNIFORM RESIDUAL DIS- 

C TRIBUTIGN (WHEN RSW=1) 

C ISEED - FI^ST OF TWO INITIAL SEEDS FOR THE RANDOM 

C NUMBER GENERATORS 

C JSfcED - SECOND OF TWO INITIAL SEEDS FOi" THE "ANDOM 

C NUMBER GENERATORS 

C 

c note: 
c 

C * THE NUMBER OF ACTIVE CONSTRAINTS IS EQUAL TO M+l t THE 

C NUMBER OF INACTIVE CONSTRAINTS IS EQUAL TO 2**(M+K+1) 

C For M «GE» 1 AND K oGEo 

C 

C * THE DEFAULT VALUE FOR EACH OF THE SWITCHES BSfl.CSW.ISW. 

C JSWtMSW»FSW IS ZE30 

c 

C * WHEN EPS=0« NO PERTURBATION OF EQUALLY SPACED OBSEF.VA- 

C TION POINTS IS PERFORMED (MSW=0). OTHERWISE, A REASONABLE 

C ~ANGE IS < EPS < 0.5 

C 

C output arguments: 

C 

C BETA SOLUTION VECTOR TO THE GENERATED PROBLEM OF 

C DIMENSION Ml 

C PR3MAT- THE GENERATED PROBLEM MATRIX WITH KKM1 ROWS 

C AND M2 COLUMNS 

C SUM3ES- SUM OF THE ABSOLUTE VALUES OF THE RESIDUALS IN THE 

C OPTIMAL LI SOLUTION 

C KKM1 — NUMBER OF .10WS IN PRBMAT 

C M2 NUMBER OF COLUMNS IN PRBMAT 

C LCC VECTOR OF DIMENSION Ml CONTAINING THE F OW LOCATIONS 

C OF THE ACTIVE CONSTRAINTS IN PROBLEM MATFIX PF.BMAT 

C ISEED - FIRST OF TWO RANDOM SEEDS AVAILABLE UPON RETURN 

C JSEED - SECOND OF TWO RANDOM SEEDS AVAILABLE UPON RETURN 
C 

c restrictions: 

c 

C M MUST BE .GE. I 

C K MUST BE »GE. 

C M+K MUST BE .LE. 8 

C NRCW MUST BE »GEo M+l +2** ( M+K+ 1 ) 

C 

C***** ************************************** *************************** 

C 

C SUBPROGRAMS USED: 

c 

C ABS FG3T3AN LIB^AnY FUNCTION 

C UNI RAN -- SUBROUTINE TO GENERATE PSEUDO-RANDOM NUMBERS FROM 
C A UNIFCi-.M DISTRIBUTION OVE? [Oil ] 

C NOF.RAN SUBROUTINE TO GENERATE PSEUDO-RANDOM NUMBERS FROM 

C A NORMAL DISTRIBUTION WITH MEAN ZERO AND STANDARD 

C DEVIATION EQUAL TO ONE 

C SORTM SUBROUTINE USED TO SORT THE ELEMENTS OF A VECTOR 

C INTO INCREASING ORDER, ALSO PROVIDES A MAPPING 

C (N=>CS) FSGM OLD POSITION TO NEW (SORTED) POSITION 
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c 


SET1 SUBROUTINE USED TO CALCULATE EQUIDISTANT OBSER- 


c 


VATION POINTS FO^ THE INACTIVE CONSTRAINTS 


c 


SET 2 SUBROUTINE USED TO CALCULATE NON-EQUIDISTANT 


c 


OBSERVATION POINTS F OF: THE INACTIVE CONSTRAINTS 


c 


EQUAL SUBROUTINE TO PRODUCE OBSERVATION POINTS EQUALLY 


c 


SPACED ON £0*1 ], THOUGH FOSSIBLY PERTURBED BY A 


c 
c 
c 
c 
c 
c 
c 


RANDOM AMOUNT IF EPS > 


LANGUAGE: ANSI FORTRAN CI 966) 


references: 


K»L«HOFFMAN & D.Re SH I E~ • • A TEST PROBLEM GENERATOR FO? 


c 


DISCRETE LINEAR LI APPROXIMATION PROBLEMS^, TO APPEAR IN 


c 


TFiANS* ON MATH» SOFTWARE 


c 


P«D«DOMICH & D.R.SHIE3* *GENE^ATO^S FOR DISCRETE POLYNOMIAL 


c 
c 
c 
c 
c 


LI APPROXIMATION PRGBLEMS i 


WRITTEN BY: 


PAUL DOMICH & DOUGLAS SHIER 


c 


CENTEf FO& APPLIED MATHEMATICS 


c 


NATIONAL BUREAU OF STANDARDS 


c 
c 
c 
c 
c* 


WASHINGTON D»C© 20234 


JANUARY, 1979 


*************************** ****^ 


c 
c 


* START OP SUBROUTINE POLY2 * 



INTEGER 8SW,CSW.RSW,LCC( 1 ) 

$EAL PRDMAT(NRGW, 1 ) ,BETA( i ) 

DIMENSION X<52 1) .A(256,2) ,UU( 521 •10)«NPGSC 521 > 

DOUBLE PRECISION HOLD • HOLD 1 . XPT , X*T 1 » k I tR2 ♦ DX * DX1 • DSUM 

DATA BLFT.BRIT/-10 .0 • 10.0/ 

DATA CLFT.CRIT/-1.0, 1.0/ 
C 

c initialization: THE GENERATED problem WILL HAVE Ml ACTIVE 
C CONSTRAINTS* KK INACTIVE CONSTRAINTS WITH POSITIVE RESIDUALS, 
C AND KK INACTIVE CONSTRAINTS WITH NEGATIVE RESIDUALS 
C 

I0UT=6 

DSUM=0.0D0 

M1=M+1 

M2=M+2 

MK=M+K 

KK-2**MK 

KK2=KK+KK 

KM1=KK+M1 

KKM1=KK2+M1 
C 

C CHECK THE CONSISTENCY OF THE IN°UT PA^AMETE^S 
C 

IF(M »LT« 1 oG^o K oLTo 0) GOTO 9990 

IFCMK *GT. 9) GOTO 9991 

IFCKKMl .GTa NROW) GOTO 9991 

IPCBLEFT .GT. BRIGHT) GOTO 1001 
1000 IF(CLEFT # GT« C^IGHT) GOTO 1002 
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GOTO 100 3 
100 1 XHOLD^BLEF T 

BLEFT=P" IGHT 
BR I GHT-XHUL13 
GOTO 1000 
1002 XHOLO=CLEFT 
CLEFT=C^IGHT 
CRIGHT=XHQLD 



1003 

C 

C 


CONTINUE 


DETERMINE INTERVALS USED 


c 
c 


VALUES 


IF(BSW .NE. 1)3LEFT=BLFT 




IF(BSW .NE. 1 ) BRIGHT^BRIT 




IF(CSW oNEo 1)CLEFT=CLFT 




IFCCS* .NE* 1 )CRIGHT=CRI T 



IN PAQBLEM By SPECIFICATION OR DEFAULT 



BINTEf =B3IGHT-BLEF T 

C INTER=CRIGHT-CLEFT 
C 

C DETERMINE THE SOLUTION VECTCR BY USER SPECIFICATION OR 3Y 

C FANDGM GENERATION 

C 

IF(JSW oNEo 1) GOTO 20 1 

CALL UNIRANCM1 9 l«Xt ISEEOt JSECOJ 

DO 2000 1=1. Ml 

BETA( I )=BLEFT+X( I )*8 INTER 

2000 CONTINUE 

2001 CONTINUE 
C 

C CALCULATE THE Ml OBSERVATION POINTS CORRESPON D I NG TO THE 
C ACTIVE CONSTRAINTS© NGTl: ACTIVE CONSTRAINTS A? E DEFINED 
C TO HAVE IDENTICALLY ZcRO RESIDUALS 
C 

IF<MSW .EQ« 1KALL UN I RAN ( M I ♦ 1 . X , I SEE J , J SE EG > 

IF(MSW oNEo 1)CALL EQU AL ( M 1 • X , EPS . I SEE D . J SEED ) 

DO 300 1 1=1 t Ml 

HOLO=BETA(M1 ) 

UU( I «2) = 1« J 

DX=lc0D0 

XPT=CLFFT+X< I )*CINTER 

DO 300 J=ltM 

N = M 1 - J 

H0LD=3ETA(N)+XPT*H0LD 

DX=XPT*DX 

UUC I »J+2)=DX 
3)00 CONTINUE 

UUC I • 1 )=HOLD 
3001 CONTINUE 
C 

C CALCULATE THE OBSERVATION POINTS CORRESPONDING TO THE INACTIVE 
C CONSTRAINTS BY "SETl* (EQUIDISTANT) OR •SET2" (NON-EQUIDISTANT) 
C 

IF(ISW .EG* 1>CALL SETl(KK.A) 

IF(ISW NE* 1)CALL SET 2 ( KK • MK . A» I SEED • USE ED ) 
C 

C CALL RANDOM NUMBER GENERATOR FOP CALCULATION OF kESIDUALS 
C AND GENLRATE THE DESIGN MATRIX FOR THE INACTIVE CONSTRAINTS 
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IFiHSU oEQ» l)CALL UNI R AN ( KK2 , 1 * X * I SEED , JSEEO ) 

IF<RSW .NE. 1>CALL NGRRAN{ KK2 , 1 ♦ X , I SEED, J SEED ) 

00 4001 I=M2,KM1 

N=I-M1 

IKK=I+KK 

NKK=N+KK 

H0LD = 8ETA(M1 ) 

H0LD1=HCL0 

UU( h2)-U0 

UU< IKK, 21=1.0 

DX=1.0D0 

DX1=1 # 0D0 

XPT=CLEFT+A(N» 1 >*C INTER 

XPT1=CLEFT+A(N,2)*CINTEF. 

DO 4000 J=ltM 

L=M1-J 

HGLO=BETA(L)+XPT*HCLD 

H0LD1=BETA(L) +XPT 1*H0LD1 

DX=XPT*DX 

DX1=XPT1*DX1 

UU( I • J+2)=DX 

UU< IKK, J+2 )=DX1 

4000 CONTINUE 
C 

C CALCULATE THE RESIDUALS AND ADD TO Y-VALUES 
C 

R1=ABS(RESID*X(N) ) 

R2 = ABS(P.ESID*X(NKK ) ) 

H0LD=HCLD+R1 

H0L01=HCLD1-?12 

DSUM = DSU>1 + R1+R2 

UU( I , 1 >=H0LD 

UU( IKK, 1) =H0LD1 

4001 CONTINUE 
SUMRES=DSUM 

C 

C SORT THE OBSERVATION POINTS, REARRANGE UU AND STORE THE RESULT 

C IN PR8MAT 

C 

DO 5000 I=l»KKMl 

X ( I )=UU( I ,3) 
5000 CONTINUE 

CALL SCRTM<X,KKM1 , NPOS ) 

DO 5002 I-1.KKM1 

KLGC=NFCS< I ) 

DO 5001 J=1«M2 

P^BMATCKLOC, J >=UU< I, J) 
500 1 CONTINUE 
5002 CONTINUE 

DO 5CC3 1=1, Ml 

LOC( I l=NPOS< I ) 
5CC3 CONTINUE 

GOTO 5999 

9990 WRITE( IOUT. 1 ) M.K 
GOTO 9999 

9991 WRITE( ICUT*2) M»K 

1 FCRMATC 10X»15HINPUT VALUE OF M = . 1 5 , 2 X. 7H0R K =,I5.2X,3HIS , 
+ 7HILLEGAL) 

2 FORMAT! 10X,20HINPUT VALUES OF M = , I 5 • 2X, 8HAND K = , 15 , 2X , 6HEXCECD 
+ ,18H AVAILABLE STORAGE) 

9959 FETUSN 
END 
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SUBROUTINE SETMKK.A) 

C 

c* ********** ********************************************************** 

c 

C THIS SUBROUTINE CREATES A SET OF 2*KK REAL NUMBERS EQUALLY 

C SPACED ON £0.13. THIS SET IS DIVIDED INTO TWC SUBSETS OF KK 

C ELEMENTS WHICH HAVE THE FOLLOWING PROPERTIES: THE SUM OF THE 

C P-TH POWERS OF ELEMENTS IN THE FIRST SUBSET EQUALS TH> SUM OF 

C THE P-TH POWERS OF ELEMENTS IN THE SECOND SUBSET. FOR VALUES 

C P = 0. |, .... L0G2CKK). 

C 

C THE VALUES OF THE ELEMENTS RANGE BETWEEN AND 1, WITH THE 

C SMALLEST VALUE IDENTICALLY ZERO AND THE LARGEST VALUE IDEN- 

C TICALLY ONE. 

C 

C THE FIRST SUBSET IS RETURNED IN THE FIRST COLUMN OF ARRAY A 

C AND THE SECOND SUBSET IS RETURNED IN THE SECOND COLUMN. 

C 

c****** * * ** *** * ** *** ****************** ***** **************************** 

c 

DIMENSION AC256.2) 

DOUBLE PRECISION XD 
C 

C***** INITIALIZE THE FIRST TWO ELEMENTS OF EACH SUBSET 
C 

A (1,1 I -0.0 

A < I . 2 > =1 . O 

A{2. 11=3.0 

A(2.2)=2.0 

N = 2 
C 

C***** CALCULATE THE REMAINING ELEMENTS OF THE TWO SUBSETS 
C 

IF<KK,EQ.2> GO TO 120 

DO 100 1=3. KK 

M=I-(2**<N-1I ) 

A( I ,1 )=A(M,2)*FL0AT(2**N) 

A(I . 2)=A<M.1)+FL0AT<2**N) 

IF <<2**(N-1)> .SO. Mi N=N+1 
100 CONTINUE 

C 

C***** STANDARDIZE THE VALUES TO LIE IN C0«1] 
C 

XD=2*KK-1 

DO 110 I=1.KK 

Ad.l J=A(I.1I/XD 

A(I,2)=A( I .2)7X0 
110 CONTINUE 
120 RETURN 

END 
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SUBROUTINE SET2 i KK.MK. A, ISEED. JSEED) 
C 

c ************************************************************* ******** 
c 

C THIS SUBROUTINE CREATES A SET OF 2*KK REAL NUMBERS UNEQUALLY 

C SPACED ON £0,13. THIS SET IS DIVIDED INTO TWO SUBSETS OF KK 

C ELEMENTS WHICH HAVE THE FOLLOWING PROPERTIES: THE SUM OF THE 

C P-TH POWERS OF ELEMENTS IN THE FIRST SUBSET EQUALS THE SUM OF 

C THE P-TH POWERS OF ELEMENTS IN THE SECOND SUBSET, FOR VALUES 

C P=0. !•*••. L0G2CKK). 

C 

C THE VALUES OF THE ELEMENTS RANGE BETWEEN AND 1, WITH THE 

C SMALLEST VALUE IDENTICALLY ZERO AND THE LARGEST VALUE IDEN- 

C TICALLY ONE. 

C 

C THE FIRST SUBSET IS RETURNED IN THE FIRST COLUMN OF ARRAY A 

C AND THE SECOND SUBSET IS RETURNED IN THE SECOND COLUMN. 

C 

C***** ****************** *********************************************** 

C 

DIMENSION A( 256,2) ,X<9> 

DOUBLE PRECISION XHIGH 

A(1.1)=0»0 

MK1=MK+1 

M = l 

1 = 1 
C 

C***** SELECT ELEMENT A(1.2> RANDOMLY 
C 

CALL UNI RAN(MK 1.1. X.I SEED. JSEEDI 

A CI. 2>= X(I) 

XHIGH=X( I ) 

1*1+1 
C 

C***** CALCULATE THE REMAINING ELEMENTS OF THE TWO SUBSETS 
C 

DO 100 J=2.KK 

A< J.l J=A(M,2I+X(I ) 

A( J.2)=ACM.l J+X(i) 

IF {XHIGH .LT. A<J.2II XHIGH=A<J.21 

IF (XHIGH .LT« ACJ.ll) XHIGH=A4J.l> 

IF (FLOAT! JI/2«0-FL0AT(MI .GE. 0*1) GOTO 90 

M = 

1 = 1+1 
90 CONTINUE 

M = M+1 
100 CONTINUE 

C 

C***** STANDARDIZE THE VALUES TO LIE IN CO.l] 
C 

DO 110 1=1. KK 

All .1 )=ACI.l l/XHIGH 

A( I,2) = A(I .2 I /XHIGH 
110 CONTINUE 

RETURN 

END 
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SUBROUTINE NORRANCN, ISTART. X. I SEED, JSEED) 
C 

C PURPOSE THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N 

C FROM THE THE NORMAL (GAUSSIAN) 

C DISTRIBUTION WITH MEAN = AND STANDARD DEVIATION = 1< 

C THIS DISTRIBUTION IS DEFINED FOR ALL X AND HAS 

C THE PROBABILITY DENSITY FUNCTION 

C FIX) = (l/SQRT(2*PIl>*EXP(-X*X/2> . 

C INPUT ARGUMENTS — N « THE DESIRED INTEGER NUMBER 

C OF RANDOM NUMBERS TO BE 

C GENERATED. 

C ISTART = AN INTEGER FLAG CODE WHICH 

C (IF SET TO OJ WILL START THE 

C GENERATOR OVER AND HENCE 

C PRODUCE THE SAME RANDOM SAMPLE 

C OVER AND OVER AGAIN 

C UPON SUCCESSIVE CALLS TO 

C THIS SUBROUTINE WITHIN A RUN; OR 

C (IF SET TO SOME INTEGER 

C VALUE NOT EQUAL TO 0. 

C LIKE, SAY, 1) WILL ALLOW 

C THE GENERATOR TO CONTINUE 

C FROM WHERE IT STOPPED 

C AND HENCE PRODUCE DIFFERENT 

C RANDOM SAMPLES UPON 

C SUCCESSIVE CALLS TO 

C THIS SUBROUTINE WITHIN A RUN. 

C OUTPUT ARGUMENTS X = A SINGLE PRECISION VECTOR 

C (OF DIMENSION AT LEAST N) 

C INTO WHICH THE GENERATED 

C RANDOM SAMPLE WILL BE PLACED. 

C OUTPUT A RANDOM SAMPLE OF SIZE N 

C FROM THE NORMAL DISTRIBUTION 

C WITH MEAN = AND STANDARD DEVIATION =1. 

C PRINTING NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 

C RESTRICTIONS THERE IS NO RESTRICTION ON THE MAXIMUM VALUE 

C OF N FOR THIS SUBROUTINE, 

C OTHER DATAPAC SUBROUTINES NEEDED UNIRAN. 

C FORTRAN LIBRARY SUBROUTINES NEEDED ALOG, SORT, SIN, COS. 

C MODE OF INTERNAL OPERATIONS SINGLE PRECISION. 

C LANGUAGE — ANSI FORTRAN. 

C METHOD — BOX-MULLER ALGORITHM. 

C REFERENCES — BOX AND MULLF^ »A NOTE ON THE GENERATION 

C OF RANDOM f urtMAL DEVIATES 1 . JOURNAL OF THE 

C ASSOCIATIOIS FOR COMPUTING MACHINERY. 1958, 

C PAGES 610-611. 

C — TOCHER. THE ART OF SIMULATION, 

C 1963, PAGES 33-34. 

C — HAMMERSLEY AND HANDSCOMB, MONTE CARLO METHODS, 

C 1964, PAGE 39. 

C — JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE 

C DISTRIBUTIONS — 1, 1970, PAGES 40-111. 

C WRITTEN BY — JAMES J. FILLIBEN 

C STATISTICAL ENGINEERING LABORATORY (205.031 

C NATIONAL BUREAU OF STANDARDS 

C WASHINGTON, D. C. 20234 

C PHONE: 30 1-921-2315 
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C ORIGINAL VERSION — JUNE 1972. 

C UPDATED — SEPTEMBER 1975. 

C UPDATED —NOVEMBER 1975. 

C UPDATED — JULY 1976. 

C 

C 

DIMENSION X(l) 

DIMENSION Y<2) 

DATA PI/3.14159265358979/ 
C 

IPR=6 
C 

C CHECK THE INPUT ARGUMENTS FOR ERRORS 
C 

IF(N.LT.lIG0TO50 

G0T090 
50 KRITEIIPR. 51 

WRITE(IPR.47)N 

RETURN 
9 CONTINUE 
5 FORMAT! I H • 91 H***** FATAL ERROR — THE FIRST INPUT ARGUMENT TO THE 
1 NORRAN SUBROUTINE IS NON-POSITIVE *****) 
47 FORMATUH • 35H***** THE VALUE OF THE ARGUMENT IS .18 .6H *****) 

C 

C START POINT 

C 

C GENERATE N UNIFORM 10.11 RANDOM NUMBERS; 

C THEN GENERATE 2 ADDITIONAL UNIFORM (0.1) RANDOM NUMBERS 

C (TO BE USED BELOW IN FORMING THE N-TH NORMAL 

C RANDOM NUMBER *H£N THE DESIRED SAMPLE SIZE N 

C HAPPENS TO QE ODD). 

C 

CALL UNIRANtN.ISTART.X.ISEED. JSEED) 

CALL UNIRAN<2, l.Y.I SEED, JSEED) 
C 

C GENERATE N NORMAL RANDOM NUMBERS 

C USING THE BOX-MULLER METHOD. 

C 



D0200I=1.N.2 

IPl=i*l 

U1=X<I) 

IF( I.EQ.N)GOT0210 

U2=X(IP1 ) 

GOTO220 
210 U2=Y(2) 
220 ARGI=-2.0*AL0GCU1> 

ARG2=2.0*PI*U2 

SQRT1-SQRT<ARG1 ) 

Z1=SQRT1*C0S( ARG2) 

Z2=SQRT1*SIN(ARG2) 

X(I l-Zl 

IF( I.EQ.N)G0TO200 

X(IP1)=Z2 
200 CONTINUE 

RETURN 
END 
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SUBROUTINE EQUAL ( M 1 . X. EPSt ISEEDt JSEEOI 

C 

c *********** ************************************************************ 

c 

C THIS SUBROUTINE GENERATES EQUALLY SPACED OBSERVATION POINTS 

C THOUGH POSSIBLY PERTURBED BY A RANDOM AMOUNT <IF SPS .NE. 0.0). 

C THE VALUES OF THESE OBSERVATION POINTS RANGE BETWEEN AND 

C It WITH THE SMALLEST VALUE IDENTICALLY ZERO AND THE LARGEST 

C VALUE IDENTICALLY ONE, THE PERTURBATION APPLIED (WHEN EPS 

C .NE. 0.0) IS GENERATED FROM A UNIFORM PROBABILITY DISTRIBUTION. 

C 

c*** ******************************** ************************************ 

C 

REAL X<1) 

DOUBLE PRECISION H»EPSA 

M=M1-I 

CALL UN I RANI Ml , ltX, I SEED, JSE2D) 

X( 1 )=0.0 

XCM1 |*1«0 

IF(M .LT. 2) GOTO 25 

H=l ,0/FLOAT(M) 

EPSA=H*EPS 

DO 1000 I=2tM 

X{ I )=H*F LOAT< 1-1 )*EPSA*<2.0*X{ I)-l .0) 
100 CONTINUE 
25 RETURN 

END 



SUBROUTINE SORTM ( X .N , NPOS) 

C 

C*** ************************************************ ****************** 

c 

C THIS SUBROUTINE SORTS THE INPUT SI NGLE-PRECI SON VECTOR X WITH 

C N ELEMENTS. AND RETURNS IN X THE ELEMENTS SORTED INTO ASCENDING 

C ORDER. THE VECTOR NPOS GIVES THE MAPPING FROM THE ORIGINAL 

C POSITION IN X TO THE NEW SORTED POSITION. THAT IS, XII) WILL 

C NOW CORRESPOND TO THE NPOSf I ) - TH SMALLEST ELEMENT OF THE 

C VECTOR X. 

C 

c*** ******************************************** ********************* 

C 

DIMENSION XC1 ),NPGS(1) .XP0S(600) 
CALL SGRTPIX,N,X.XPOS) 
00 10 1=1, N 
NX=XPOS( I )*0. 00005 
NPOS(NX)-! 
10 CONTINUE 
RETURN 
END 
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SUBROUTINE SORTPI X *N. Y* XPOS1 
C 

C PURPOSE — THIS SUBROUTINE SORTS (IN ASCENDING ORDER) 

C THE N ELEMENTS OF THE SINGLE PRECISION VECTOR X, 

C PUTS THE RESULTING N SORTED VALUES INTO THE 

C SINGLE PRECISION VECTOR Y* 

C AND PUTS THE POSITION UN THE ORIGINAL VECTOR X) 

C OF EACH OF THE SORTED VALUES. 

C INTO THE SINGLE PRECISION VECTOR XPOS. 

C THIS SUBROUTINE GIVES THE DATA ANALYST 

C NOT ONLY THE ABILITY TO DETERMINE 

C WHAT THE MIN AND MAX IFOR EXAMPLE) 

C OF THE DATA SET ARE, BUT ALSO 

C WHERE IN THE ORIGINAL DATA SET 

C THE MIN AND MAX OCCUR. 

C THIS IS ESPECIALLY USEFUL FOR 

C LARGE DATA SETS. 

C INPUT ARGUMENTS — X = THE SINGLE PRECISION VECTOR OF 

C OBSERVATIONS TO BE SORTED. 

C N = THE INTEGER NUMBER OF OBSERVATIONS 

C IN THE VECTOR X. 

C OUTPUT ARGUMENTS — Y = THE SINGLE PRECISION VECTOR 

C INTO WHICH THE SORTED DATA VALUES 

C FROM X WILL &E PLACED. 

C — XPOS ■ THE SINGLE PRECISION VECTOR 

C INTO WHICH THE POSITIONS 

C UN THE ORIGINAL VECTOR X) 

C OF THE SORTED VALUES 

C WILL BE PLACED. 

C OUTPUT — THE SINGLE PRECISION VECTOR XS 

C CONTAINING THE SORTEO 

C (IN ASCENDING ORDER) VALUES 

C OF THE SINGLE PRECISION VECTOR X. AND 

C THE SINGLE PRECISION VECTOR XPC S 

C CONTAINING THE POSITIONS 

C <IN THE ORIGINAL VECTOR X) 

C OF THE SORTED VALUES. 

C PRINTING NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 

C RESTRICTIONS THE DIMENSIONS OF THE VECTORS IL AND IU 

C (DEFINED AND USED INTERNALLY WITHIN 

C THIS SUBROUTINE) DICTATE THE MAXIMUM 

C ALLOWABLE VALUE OF N FOR THIS SUBROUTINE. 

C IF IL AND IU EACH HAVE DIMENSION K, 

C THEN N MAY NOT EXCEED 2**<K+1) -I. 

C FOR THIS SUBROUTINE AS WRITTEN* THE DIMENSIONS 

C OF IL AND IU HAVE BEEN SET TO 36* 

C THUS THE MAXIMUM ALLOWABLE VALUE OF N IS 

C APPROXIMATELY 137 BILLION. 

C SINCE THIS EXCEEDS THE MAXIMUM ALLOWABLE 

C VALUE FOR AN INTEGER VARIABLE IN MANY COMPUTERS* 

C AND SINCE A SORT OF 137 BILLION ELEMENTS 

C IS PRESENTLY IMPRACTICAL AND UNLIKELY* 

C THEN THERE IS NO PRACTICAL RESTRICTION 

C ON THE MAXIMUM VALUE OF N FOR THIS SUBROUTINE. 

C CIN LIGHT OF THE ABOVE* NO CHECK OF THF 

C UPPER LIMIT OF N HAS BEEN INCORPORATED 

C INTO THIS SUBROUTINE.) 
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C OTHER DATAPAC SUBROUTINES NSEDEO — NONE. 

C FORTRAN LIBRARY SUBROUTINES NEEDED NONE. 

C MODE OF INTERNAL OPERATIONS SINGLE PRECISION, 

C LANGUAGE ANSI FORTRAN. 

C COMMENT — THE SMALLEST ELEMENT OF THE VECTOR X 

C WILL BE PLACED IN THE FIRST POSITION 

C OF THE VECTOR Y, 

C THE SECOND SMALLEST ELEMENT IN THE VICTOR X 

C WILL BE PLACED IN THE SECOND POSITION 

C OF THE VECTOR Y. 

C ETC. 

C COMMENT THE POSITION (I THROUGH N) IN X 

C OF THE SMALLEST ELEMENT IN X 

C WILL BE PLACED IN THE FIRST POSITION 

C OF THE VECTOR XPOS. 

C THE POSITION Cl THROUGH N) IN X 

C OF THE SECOND SMALLEST ELEMENT IN X 

C WILL BE PLACED IN THE SECOND POSITION 

C OF THE VECTOR XPOS. 

C ETC. 

C ALTHOUGH THESE POSITIONS ARE NECESSARILY 

C INTEGRAL VALUES FROM 1 TO N. IT IS TO B5 

C NOTED THAT THEY ARE OUTPUTED AS SINGLE 

C PRECISION INTEGERS IN THE SINGLE PRECISION 

C VECTOR XPOS. 

C XPOS IS SINGLE PRECISION SO AS TO BE 

C CONSISTENT WITH THE FACT THAT ALL 

C VECTOR ARGUMENTS IN ALL OTHER 

C DATAPAC SUBROUTINES ARE SINGLZ PRECISION. 

C COMMENT THE INPUT VECTOR X REMAINS UNALTERED. 

C COMMENT IF THE ANALYST DESIRES A SORT *IN PLACE* • 

C THIS IS DONE BY HAVING THE SAME 

C OUTPUT VECTOR AS INPUT VECTOR IN THE CALLING SEQUENCE 

C THUS* FOR EXAMPLE. THE CALLING SEQUENCE 

C CALL SORTF{X.N.X.XPOS> 

C IS ALLOWABLE AND WILL RESULT IN 

C THE DESIRED •IN-PLAC£« SORT. 

C COMMENT THE SORTING ALGORTHM USED HEREIN 

C IS THE BINARY SORT. 

C THIS ALGORTHIM IS EXTREMELY FAST AS THE 

C FOLLOWING TIME TRIALS INDICATE. 

C THESE TIME TRIALS WERE CARRIED OUT ON THZ 

C UNI VAC 1108 EXEC 8 SYSTEM AT NBS 

C IN AUGUST OF 1974. 

C BY WAY OF COMPARISON. THE TIME TRIAL VALUES 

C FOR THE EASY-TO-PROGRAM BUT EXTREMELY 

C 

c 

C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
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INEFFICIENT BUBBLE 


SORT 


ALGORITHM HAVE 






ALSO BEEN INCLUDED- 


- 










NUMBER OF RANDOM 




BINARY SOFT 






BUBBLE SO 


NUMBERS SORTED 












N s 10 




•002 SEC 






.002 SEC 


N as 100 




•Oil SEC 






.045 SEC 


N » 1000 




.141 SEC 






4.33 2 SEC 


N = 3000 




.476 SEC 






37.683 SEC 


N * 10000 




1.887 SEC 




NOT COMPUTED 


REFERENCES — CACM MARCH 1969* 


PAGE 186 (BINARY 


SORT 


ALGORITHM 


BY RICHARD C. SINGLETON). 








CACM JANUARY 1970. PAGE 54. 









C CACM OCTOBER 1970, PAGE 624. 

C JACM JANUARY 1961* PAGE 41* 

C WRITTEN BY — JAMES J. FILL I SEN 

C STATISTICAL ENGINEERING LABORATORY 1205.031 

C NATIONAL BUREAU OF STANDARDS 

C WASHINGTON, D. C. 20234 

C PHONE 30 I -92 1-23 1 5 

C ORIGINAL VERSION JUNE 1972. 

C UPDATED — NOVEMBER 1975* 

C 

C 

DIMENSION X<1).YI1>,XP0S< 1) 

DIMENSION IU(36),IL(36) 
C 

C CHECK THE INPUT ARGUMENTS FOR ERRORS 

C 

IPR=6 

IFCN.LT.1IG0T050 

IF(N.EQ.1)G0T055 

H0LD=X<1) 

D060I=2,N 

I F < X < I ) . NE . HOLD ) GOTO 90 

60 CONTINUE 
WRITE<IPR» 9)H0LD 
D061 1=1 ,N 
YCII=X<II 
XPOSCI )=I 

61 CONTINUE 
RETURN 

50 WRITEIIPR,15) 
WRITE<IPR.47)N 
RETURN 
55 WRITE! IPR,18I 
Y( 1I=X(1) 
XPOSd ) = 1.0 
RETURN 
90 CONTINUE 
9 FORMAT! 1H ,108H***** NON-FATAL DIAGNOSTIC—THE FIRST INPUT ARGUME 
INT I A VECTOR) TO THE SORTP SUBROUTINE HAS ALL ELEMENTS « .El 5.8. 6 
1H *****) 

15 FORMAT! IH . 91H***** FATAL ERROR THE SECOND INPUT ARGUMENT TO THE 

1 SORTP SUBROUTINE IS NON-POSITIVE *****> 

18 FORMATCIH , 100H***** NON-FATAL DIAGNOSTIC THE SECOND INPUT ARGUME 

INT TO THE SORTP SUBROUTINE HAS THE VALUE 1 *****) 
47 FORMATCIH • 35H***** THE VALUE OF THE ARGUMENT IS ,18 • 6H *****) 

C 

C START POINT 

C 

C COPY THE VECTOR X INTO THE VECTOR Y 
DOlOOI^l.N 
YC I) = X< I) 
100 CONTINUE 
C 

C DEFINE THE XPOS {POSITION) VECTOR. BEFORE SORTING, THIS WILL 
C BE A VECTOR WHOSE I-TH ELEMENT IS EQUAL TO I. 
C 



D0150I=1,N 
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XPOSI I)-I 
150 CONTINUE 
C 

C CHECK TO SEE IF THE INPUT VECTOR IS ALREADY SORTED 
C 

NM1=N-1 

D0200I=1,NM1 

IP1 = H-1 

IF{YC I).LE.Y( IP1 ) JG0T0200 

G0T0250 
200 CONTINUE 

RETURN 
250 M=l 

1=1 

J=N 
305 IF(I.GE*J>GOTO370 
310 K=I 

MID=C I*J>/2 

AMED=YCMID ) 

BMED^XPOSIMID) 

IFCYC ll.L£«AMZD)G0TQ320 

Y(MID*=YCI1 

XP0SCMID)=XP0S( I ) 

Y< I >=AMED 

XPOSC I )=BMED 

AMED=YCMID) 

BMED=XPOS(MIDI 
320 L=J 

IFCYC J>«GE.AMED)GOT0 340 

YCMIDI^YC J) 

XPOSC MID >=XPOS( J) 

YC J)=AMED 

XPOSC J) = BMED 

AMED=Y(MID) 

BMED=XPOS(MIDI 

IF(YC I).LE«AMED>GOT0340 

YCMID )=Y( I ) 

XPOSCMID >=XPOS< I ) 

YCI )=AMED 

XPOSC I ) = BMED 

AMED=Y(MID> 

BMED^XPOSCMID) 

GOTO340 
330 Y(L)=YCK> 

XPOS(i.l=XPOSCK) 

YCK)=TT 

XP0SCK)=ITT 
340 L=L-1 

IF<YCt-J.GT.AMED)GOT0340 

TT=YCLI 

ITT=XP0SCL) 
350 K = K + 1 

IF(Y(KKLT.AMED)G0T0350 

IF(K,LE«L)G0T0330 

LMI=L-I 

JMK=J-K 

IF(LMI •LE.JMK)GOT0 360 

IL(M)=I 
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I=K 

M=M+1 
GOT0380 
360 ILCM|=K 
IUCM)=J 

J=L 

M=M+1 
G0T0380 
370 M = M-1 

IFCM.EQ.OJRETURN 
I=IL<M) 

380 JMI=J-I 

IF< JMI«GE»11)G0T0310 

IF<I*gQ«l)GOT0305 

1=1-1 
390 1=1+1 

IFCI.EQ.J)G0TO370 

AMEDsYCI+1 ) 

BMED=XPOS< 1 + 11 

IFiYi II.LS»AMED)GQT0390 

K=I 
395 Y<K+1)=Y(KJ 

XPOSCK+1 )=XP0SIKJ 

K=K-1 

IF<AMED,LT«Y<K> IGOTQ395 

YCK+1 )=AMED 

XPOSCK+1 )=BMEO 

GOT 03 90 

END 
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